function f=NUFFT2D(alpha,omega,precision,M)

%
% Yoel Shkolnisky, January 2010.

[b,m,q]=NUFFTParameters(precision);

if length(size(omega))~=2
    error('omega must be a 2 dimensional array')
end

if size(omega,2)~=2
    error('omega must be an array with 2 columns');
end
mM=m*M;
alpha=alpha(:);
n=numel(alpha);


if size(omega,1)~=n
    error('alpha and omega must have the same number of rows');
end

%N=4096;
%save_memory=0;
%if (2*N+1)*(N+1) < n
%	save_memory=1;
%end

%if ~save_memory
	mu=round(omega.*m);
	% All precomputations are borken into 1D ones.
	Px=zeros(q+1,n,precision);
	Py=zeros(q+1,n,precision);
	for j=-q/2:q/2
		tmp1=(m.*omega(:,1)-(mu(:,1)+j)).^2;
		tmp2=(m.*omega(:,2)-(mu(:,2)+j)).^2;
		Px(j+q/2+1,:)=exp(-tmp1/(4*b))./(2*sqrt(b*pi));
		Py(j+q/2+1,:)=exp(-tmp2/(4*b))./(2*sqrt(b*pi));
	end
	clear tmp1 tmp2
	threadN=feature('numCores');
	if strcmp(precision,'single')
		threadN=2;
		tau=NUFFT2DauxiliaryMEXMultiThreadedSinglePrecision(threadN,n,M,m,q,mu,Px,Py,alpha);
		%tau=NUFFT2DauxiliaryMEXSingleThreadedSinglePrecision(n,M,m,q,mu,Px,Py,alpha);
	elseif strcmp(precision,'double')
		tau=NUFFT2DauxiliaryMEXMultiThreadedDoublePrecision(threadN,n,M,m,q,mu,Px,Py,alpha);
		%tau=NUFFT2DauxiliaryMEXSingleThreadedDoublePrecision(n,M,m,q,mu,Px,Py,alpha);
	else
		error('precision must be single or double');
	end
	clear Px Py
%else
	% n is very big, call mex that computes Px,Py on the fly thus saves the memory they require
%	tau=NUFFT2DauxiliaryMEXSingleThreadedSinglePrecisionMemorySave(n,M,m,q,omega,b,alpha);
%end

T=fftshift(ifftn(ifftshift(tau)));
T=T*mM*mM;

low_idx_M=-ceil((M-1)/2);
high_idx_M=floor((M-1)/2);
idx=low_idx_M:high_idx_M;
E=exp(b*(2*pi*idx./(mM)).^2);
offset1=ceil((mM-1)/2);
offset2=offset1+low_idx_M+1;
E2=reshape(E'*E,[M M]);
f=T(offset2:offset2+M-1,offset2:offset2+M-1).*E2;
